program hop_exct
	implicit NONE
	
	integer(4), parameter		:: PP = 100
	integer(4), parameter		:: LM = 4
	integer(4), parameter		:: LX = 3
	integer(4), parameter		:: LD = LX*LX
	
	integer(4), parameter		:: XYMID(2) = [LX/2+1,LX/2+1]

	real(8)						:: t
	real(8)						:: tau
	real(8)						:: UU
	
	integer(4)					:: upbound(LX)
	integer(4)					:: dnbound(LX)
	
	real(8)						:: Hm(LD,LD)
	real(8)						:: evals(LD)
	
	integer(4)					:: i, j
	integer(4)					:: ix, iy
	integer(4)					:: jx, jy
	
	integer(4)					:: lwork = 3*LD
	real(8), allocatable		:: wrk(:)
	integer(4)					:: sts
	
	real(8)						:: Z(LD)
	real(8)						:: ss(LD)
	real(8)						:: sv(LX,LX)
	
	real(8)						:: emin, zz
	
	integer(4)					:: iunit, ip
	
	
!!===========================================================================
	UU =  -1.0d0
	tau = 1.0d0	
	
	forall( ix = 1:LX-1) upbound(ix) = ix + 1
	upbound(LX) = 1
	forall( ix = 2:LX) dnbound(ix) = ix - 1
	dnbound(1) = LX
	
	Hm = 0.0d0
	
	do i = 1, LD
		do j = 1, LD
			call ToCoords(i,ix,iy)
			call ToCoords(j,jx,jy)
			
			if (( i == j).and.(ix==XYMID(1)).and.(iy==XYMID(2))) Hm(i,j) = UU
			
			if (( ix == jx).and.( iy == upbound(jy))) Hm(i,j) = tau
			if (( ix == jx).and.( iy == dnbound(jy))) Hm(i,j) = tau
			
			if (( iy == jy).and.( ix == upbound(jx))) Hm(i,j) = tau
			if (( iy == jy).and.( ix == dnbound(jx))) Hm(i,j) = tau
		end do
	end do
	
	allocate( wrk(lwork), stat = sts)
	call DSYEV( 'V', 'U', LD, Hm, LD, evals, wrk, lwork, sts)
	deallocate( wrk)
	
		t = 1.0d0
		Z(:) = exp(- evals(:) / t)
		zz = sum( Z(:))
		print *, LM, zz
		
		forall( i = 1:LD) Z(i) = expserL(- evals(i) / t, LM)
		zz = sum( Z(:))
		print *, zz
	stop

	emin = minval( evals)
	evals = evals - emin
	
	open( newunit = iunit, file = 'szL.txt', action = 'write', status = 'replace')
	
	do ip = 1, PP
		
		t = 0.1d0 * ip
		
		!Z(:) = exp(- evals(:) / t)
		forall( i = 1:LD) Z(i) = expserL(- evals(i) / t, LM)
		zz = sum( Z(:))
	
		sv = 0.0d0
		do i = 1, LD
			call ToCoords(i,ix,iy)
			sv(ix,iy) = sum( Z(:) * Hm(i,:)**2)
		end do
	
		sv = sv / zz
		
		write( iunit, '(F8.2,2F10.4)') t, sv(XYMID(1),XYMID(2)), sv(upbound(XYMID(1)),XYMID(2))
		
	end do
	
	close( iunit)
	
	open( newunit = iunit, file = 'sz.txt', action = 'write', status = 'replace')
	
	do ip = 1, PP
		
		t = 0.1d0 * ip
		
		Z(:) = exp(- evals(:) / t)
		zz = sum( Z(:))
	
		sv = 0.0d0
		do i = 1, LD
			call ToCoords(i,ix,iy)
			sv(ix,iy) = sum( Z(:) * Hm(i,:)**2)
		end do
	
		sv = sv / zz
		
		write( iunit, '(F8.2,2F10.4)') t, sv(XYMID(1),XYMID(2)), sv(upbound(XYMID(1)),XYMID(2))
		
	end do
	
	close( iunit)
	
	stop
	
contains
	
pure subroutine ToCoords( n, ix, iy)
	! Variables
	integer(4), intent(IN)			:: n
	integer(4), intent(OUT)			:: ix, iy

	! Body of ToCoords

	ix = ( n - 1) / LX + 1
	iy = mod( n - 1, LX) + 1

end subroutine ToCoords

pure function expSerL( x, LL)
	! Variables
	real(8)					:: expSerL
	real(8), intent(IN)		:: x
	integer(4), intent(IN)	:: LL
	!----------------------------------------------
	integer(4)			:: ip
	real(8)				:: y, z

	! Body of expSerL
	
	y = 1.0d0
	z = 1.0d0
	do ip = 1, LL
		z = z * x / ip
		y = y + z
	end do
	expSerL = y	

end function expSerL

end program hop_exct